%     AERO LOAD - программа для нагружения агрегатов самолета распределенной нагрузкой
%     Copyright (C) 2018-2024  Михаил Козлов
%     Сайт проекта: <http://www.fastenerdoctor.ru/>
%     Адрес электронной почты для связи с разработчиком:
%     <doctor.fastener@gmail.com>
% 
%     Эта программа является свободным программным обеспечением: Вы можете
%     распространять ее и/или изменять, соблюдая условия Универсальной общественной
%     лицензии GNU, опубликованной Фондом свободного программного обеспечения,
%     либо редакции 3, либо (на Ваше усмотрение) любой редакции, выпущенной позже.
% 
%     Эта программа распространяется в надежде на то, что она окажется полезной,
%     но БЕЗ КАКИХ-ЛИБО ГАРАНТИЙ, включая даже подразумеваемую гарантию ТОВАРНОЙ ПРИГОДНОСТИ,
%     либо ПРИГОДНОСТИ ДЛЯ ОПРЕДЕЛЕННЫХ ЦЕЛЕЙ. Ознакомьтесь с Универсальной общественной
%     лицензией GNU для получения более подробной информации.
% 
%     Вы должны были получить копию Универсальной общественной лицензии GNU
%     вместе с этой программой. Если Вы ее не получили, то перейдите по адресу:
%     <http://www.gnu.org/licenses/>.
% 
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <https://www.gnu.org/licenses/>.

function [CID_DISTR_LOAD_FIELD, PRESSURE_FIELD, Load] = Loading_Process(Load, Need_Reading_BDF_File, Settings_Source)

global ErrCode;

global DB_ELEMENTS;
global DB_NODES;

global ELEMENT_AREAS;
global PROJECTED_NODES;

global Cross_Product_Vector_zero_tolerance;

global ELEMENT_FORCES;
global ELEMENT_NORMAL_PROJECTIONS;

CID_DISTR_LOAD_FIELD = [];
PRESSURE_FIELD = [];

% Проверка памяти

if (strcmp(Need_Reading_BDF_File, 'FALSE'))
    if (isempty(DB_ELEMENTS))
        fprintf (1, 'В памяти программы отсутствует загруженная модель! Измените параметр Need_Reading_BDF_File.\n');
        ErrCode = -1;
        return;
    end
end

% Источник настроек нагружения

switch Settings_Source
    case 'PROGRAM_TEXT'
        fprintf (1, 'Источник настроек нагружения:                             ТЕКСТ ПРОГРАММЫ\n\n');
    case 'CONFIG_FILE'
        fprintf (1, 'Источник настроек нагружения:                             КОНФИГУРАЦИОННЫЙ ФАЙЛ\n\n');
end

% Проверка элементов

fprintf (1, 'Проверка элементов...                                     ');

for i=1:size(DB_ELEMENTS,1)
    
    switch DB_ELEMENTS{i}.type
        
        case 'Quad'
            if (~((DB_ELEMENTS{i}.N_nodes == 4) || (DB_ELEMENTS{i}.N_nodes == 8)))
                fprintf (1, '\nНедопустимое количество узлов [%i] у элемента с номером %i типа %s !\n', DB_ELEMENTS{i}.N_nodes, DB_ELEMENTS{i}.ID, DB_ELEMENTS{i}.type);
                fprintf(1, 'В программе поддерживаются Quad элементы с 4 или 8 узлами.\n');
                ErrCode = -1;
                return;
            end
            
        case 'Tria'
            if (~((DB_ELEMENTS{i}.N_nodes == 3) || (DB_ELEMENTS{i}.N_nodes == 6)))
                fprintf (1, '\nНедопустимое количество узлов [%i] у элемента с номером %i типа %s !\n', DB_ELEMENTS{i}.N_nodes, DB_ELEMENTS{i}.ID, DB_ELEMENTS{i}.type);
                fprintf(1, 'В программе поддерживаются Tria элементы с 3 или 6 узлами.\n');
                ErrCode = -1;
                return;
            end
            
        otherwise
            fprintf (1, '\nНедопустимый тип элемента с номером %i [%s] !\n', DB_ELEMENTS{i}.ID, DB_ELEMENTS{i}.type);
            fprintf(1, 'В bdf файле с моделью для нагружения должны содержаться только элементы типов Quad и Tria.\n');
            ErrCode = -1;
            return;
            
    end
    
end

fprintf (1, 'OK!\n');

% Уравнение базовой плоскости

fprintf (1, 'Получение уравнения базовой плоскости...                  ');

switch Load.Basic_Plane.Definition
    case 'POINT_NORMAL'
        Load.Basic_Plane.Equation = Plane_Equation_Point_Normal(Load.Basic_Plane.Z0, Load.Basic_Plane.Z1);
    case 'THREE_POINTS'
        Load.Basic_Plane.Equation = Plane_Equation_3_Points(Load.Basic_Plane.D1, Load.Basic_Plane.D2, Load.Basic_Plane.D3);
    otherwise
        fprintf (1, '\nОшибка функции Loading_Process: недопустимое значение параметра Basic_Plane_Definition [%s] !\n', Load.Basic_Plane.Definition);
        ErrCode = -1;
        return;
end

fprintf (1, 'OK!\n');

% Вектор действия нагрузки

Load.Vector = [Load.Basic_Plane.Equation(1); Load.Basic_Plane.Equation(2); Load.Basic_Plane.Equation(3)];
Load.Vector = Load.Vector / norm(Load.Vector);

% Проецирование угловых точек поверхности хорд на базовую плоскость

fprintf (1, 'Проецирование поверхности хорд на базовую плоскость...    ');

Load.Chord_Surface.A1_projected = Project_Point_To_Plane(Load.Chord_Surface.A1, Load.Basic_Plane.Equation);
Load.Chord_Surface.A2_projected = Project_Point_To_Plane(Load.Chord_Surface.A2, Load.Basic_Plane.Equation);
Load.Chord_Surface.A3_projected = Project_Point_To_Plane(Load.Chord_Surface.A3, Load.Basic_Plane.Equation);
Load.Chord_Surface.A4_projected = Project_Point_To_Plane(Load.Chord_Surface.A4, Load.Basic_Plane.Equation);

if (Equiv_Check_Cube(Load.Chord_Surface.A1_projected, Load.Chord_Surface.A2_projected))
    fprintf (1, '\nОшибка: проекции точек A1 и A2 на базовую плоскость совпадают!\n');
    ErrCode = -1;
    return;
end

if (Equiv_Check_Cube(Load.Chord_Surface.A1_projected, Load.Chord_Surface.A3_projected))
    fprintf (1, '\nОшибка: проекции точек A1 и A3 на базовую плоскость совпадают!\n');
    ErrCode = -1;
    return;
end

if (Equiv_Check_Cube(Load.Chord_Surface.A1_projected, Load.Chord_Surface.A4_projected))
    fprintf (1, '\nОшибка: проекции точек A1 и A4 на базовую плоскость совпадают!\n');
    ErrCode = -1;
    return;
end

if (Equiv_Check_Cube(Load.Chord_Surface.A2_projected, Load.Chord_Surface.A3_projected))
    fprintf (1, '\nОшибка: проекции точек A2 и A3 на базовую плоскость совпадают!\n');
    ErrCode = -1;
    return;
end

if (Equiv_Check_Cube(Load.Chord_Surface.A2_projected, Load.Chord_Surface.A4_projected))
    fprintf (1, '\nОшибка: проекции точек A2 и A4 на базовую плоскость совпадают!\n');
    ErrCode = -1;
    return;
end

if (Equiv_Check_Cube(Load.Chord_Surface.A3_projected, Load.Chord_Surface.A4_projected))
    fprintf (1, '\nОшибка: проекции точек A3 и A4 на базовую плоскость совпадают!\n');
    ErrCode = -1;
    return;
end

Point_1 = Load.Chord_Surface.A1_projected;
Point_2 = Load.Chord_Surface.A2_projected;
Point_3 = Load.Chord_Surface.A3_projected;

if (norm(cross((Point_2 - Point_1)/norm(Point_2 - Point_1), (Point_3 - Point_1)/norm(Point_3 - Point_1))) < Cross_Product_Vector_zero_tolerance)
    fprintf (1, '\nОшибка: проекции точек A1, A2 и A3 на базовую плоскость лежат на одной прямой!\n');
    ErrCode = -1;
    return;
end

Point_1 = Load.Chord_Surface.A2_projected;
Point_2 = Load.Chord_Surface.A3_projected;
Point_3 = Load.Chord_Surface.A4_projected;

if (norm(cross((Point_2 - Point_1)/norm(Point_2 - Point_1), (Point_3 - Point_1)/norm(Point_3 - Point_1))) < Cross_Product_Vector_zero_tolerance)
    fprintf (1, '\nОшибка: проекции точек A2, A3 и A4 на базовую плоскость лежат на одной прямой!\n');
    ErrCode = -1;
    return;
end

Point_1 = Load.Chord_Surface.A3_projected;
Point_2 = Load.Chord_Surface.A4_projected;
Point_3 = Load.Chord_Surface.A1_projected;

if (norm(cross((Point_2 - Point_1)/norm(Point_2 - Point_1), (Point_3 - Point_1)/norm(Point_3 - Point_1))) < Cross_Product_Vector_zero_tolerance)
    fprintf (1, '\nОшибка: проекции точек A1, A3 и A4 на базовую плоскость лежат на одной прямой!\n');
    ErrCode = -1;
    return;
end

Point_1 = Load.Chord_Surface.A4_projected;
Point_2 = Load.Chord_Surface.A1_projected;
Point_3 = Load.Chord_Surface.A2_projected;

if (norm(cross((Point_2 - Point_1)/norm(Point_2 - Point_1), (Point_3 - Point_1)/norm(Point_3 - Point_1))) < Cross_Product_Vector_zero_tolerance)
    fprintf (1, '\nОшибка: проекции точек A1, A2 и A4 на базовую плоскость лежат на одной прямой!\n');
    ErrCode = -1;
    return;
end

fprintf (1, 'OK!\n');

% Проецирование узлов на базовую плоскость

fprintf (1, 'Проецирование конечных элементов на базовую плоскость...  ');

PROJECTED_NODES = zeros(size(DB_NODES,1),3);

for i=1:size(DB_NODES,1)
    PROJECTED_NODES(i,:) = Project_Point_To_Plane(DB_NODES{i}.Location, Load.Basic_Plane.Equation)';
end

fprintf (1, 'OK!\n');

% Вычисление площадей элементов

if (isempty(ELEMENT_AREAS))
    
    fprintf (1, 'Вычисление площадей элементов...                          ');

    ELEMENT_AREAS = zeros(size(DB_ELEMENTS,1),1);

    for i=1:size(DB_ELEMENTS,1)

        ELEMENT_AREAS(i,1) = Get_Element_Area(i);
        
        if (ErrCode < 0)
            return;
        end

    end
    
    fprintf (1, 'OK!\n');

end

% Предварительные вычисления

fprintf (1, 'Предварительные вычисления...                             ');

Load = Precompute_Load_Data(Load);
Load = Normalize_Distributions(Load);

fprintf (1, 'OK!\n');

% Вычисление нагрузок на элементы

fprintf (1, 'Вычисление нагрузок на элементы...                        ');

ELEMENT_FORCES = zeros(size(DB_ELEMENTS,1),1);

switch Load.Distributed_Load_Mode
    
    case 'CID_DISTR_LOAD'
        
        CID_DISTR_LOAD_FIELD = zeros(size(DB_ELEMENTS,1),4);
        ELEMENT_NORMAL_PROJECTIONS = zeros(size(DB_ELEMENTS,1),1);
        
        for i=1:size(DB_ELEMENTS,1)

            ELEMENT_FORCES(i,1) = Get_Element_Force(Load, i);

            if (ErrCode < 0)
                return;
            end

            CID_DISTR_LOAD_FIELD(i,1) = DB_ELEMENTS{i}.ID;
            CID_DISTR_LOAD_FIELD(i,1+Load.Vector_Field_Component) = ELEMENT_FORCES(i,1) / ELEMENT_AREAS(i,1);
            
            Element_Normal = Get_Shell_Normal(DB_ELEMENTS{i}.ID);
            ELEMENT_NORMAL_PROJECTIONS(i,1) = dot(Element_Normal, Load.Vector);

        end
        
    case 'PRESSURE'
        
        ELEMENT_NORMAL_PROJECTIONS = zeros(size(DB_ELEMENTS,1),1);
        PRESSURE_FIELD = zeros(size(DB_ELEMENTS,1),2);
        
        for i=1:size(DB_ELEMENTS,1)

            ELEMENT_FORCES(i,1) = Get_Element_Force(Load, i);

            if (ErrCode < 0)
                return;
            end
            
            Element_Normal = Get_Shell_Normal(DB_ELEMENTS{i}.ID);
            ELEMENT_NORMAL_PROJECTIONS(i,1) = dot(Element_Normal, Load.Vector);
            
            PRESSURE_FIELD(i,1) = DB_ELEMENTS{i}.ID;
            
            if (ELEMENT_NORMAL_PROJECTIONS(i,1) == 0)
                PRESSURE_FIELD(i,2) = 0;
            else
                PRESSURE_FIELD(i,2) = - ELEMENT_FORCES(i,1) / ELEMENT_AREAS(i,1) / ELEMENT_NORMAL_PROJECTIONS(i,1);
            end
            
        end
    
end

fprintf (1, 'OK!\n');
